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Abstract 



Previous articles in this series presented a (weighted) field energy (i.e., Dirichlet integral) based 
approach to finding point source solutions to Laplace's equation over specific bounded and unbounded 
O ' domains, where the sources are assumed to be in the complimentary region. The associated mathematical 

2"^, framework was labeled a Dirichlet-integral dual-access collocation-kernel space (DIDACKS). Specifically, 

^^' for R half-space and the exterior of an M sphere, which are the primary settings used in geoexploration 

^ ' and physical geodesy, the DIDACKS approach yields exact closed-form linear equation sets for the 

,y , strengths of point sources when their locations are fixed. By building on the field energy minimization 

C/3 ' underpinnings of DIDACKS theory and by making certain natural assumptions about the general nature 

>^j^| of the energy/density configuration of the Earth's interior it is shown that the problem of estimating the 

Earth's interior density, either globally or locally, can be naturally retrained as a energy minimization 
one. Although there are certain conceptual complications to be factored in, the basic idea is that a 
static stable density configuration is a minimum energy configuration, which tends to be unique (when 
all other things are equal); hence, a field energy minimization approach can be counted on to generally 
^ ' lead to a physically motivated unique solution. Techniques touched on here should provide practical 

C^ , implementation tools, or at least some helpful hints, for handling many of the well-known ill-posedness 

issues associated with mass density estimation and other related inverse-source Laplacian problems. 
These and additional associated considerations directly lead to the possibility of flexible and powerful 
implementations of (local) mass density estimation software that incorporates and extends certain long 
^^ ' accepted and commonly used (geoexploration) implementations. Clearly, the presented techniques can 

^—J , also be directly adapted for use in other areas of applied mathematics as well as other physical problem 

areas, such as electrostatics. This article focuses more on overall implementation issues than on concrete 
specific examples and contains no numerical examples; however, due consideration has been given to the 
potential viability of the suggested approaches. 
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(i) Preamble 

This article discusses Dirichlet-Integral Based Dual- Access Collocation-Kernel (DIDACKS) based techniques 
for geoexploration or Laplacian inverse source theory. Most of the same techniques can (and have) been 
used for gravity modeling. 

From a first-cut intuitive perspective, given that physical systems in static equilibrium (or even ones in 
quasi-static equilibrium like the Earth, which, for simplicity of exposition, is treated here as a static time- 
independent case) tend to be in minimum energy configurations, and given that DIDACKS fits minimize 
the (weighted) energy of the entire external field exactly (as well as that of the exterior energy of the 
error field) and, moreover, that DIDACKS theory was first explicitly developed and tested for the two 
standard geometries that are the most commonly used ones in geoexploration and physical geodesy — namely 
the exterior of an M"^ sphere and R'^ half-space — one might tentatively conclude that in order to solve a 
geophysical inverse density estimation problem all that is necessary to do is to simply perform a DIDACKS 
point mass fit and then reinterpret the results. Here, of course, some sort of natural accommodations to 
the innate ill-conditioning of the problem must be made by, say, carefully choosing the source spacing and 
depths. On second glance, it might appear that this first take is entirely too nieve because the mathematical 
form of the gravitational field energy itself is negative (that is to say, it is proportional to the negative of 
the Dirichlet integral of the underlying potential) and thus for gravitational source estimation problems a 
negative energy solution would require minimizing the negative of the DIDACKS cost function, which would 
lead to negative run- away solutions that correspond to a worst fit rather than a best fit. It would thus seem 
that while DIDACKS theory may be directly applied as it stands to electrostatics or other problems where 
the energy density is positive and it may be applied to gravitational modeling problems (as discussed at 
length in |10j). it should not be applied to gravitational source problems without some sort of significant 
modifications. On further examination, one might be inclined to think that a direct energy minimization 
approach to gravity source estimation is more or less hopeless since mass density estimation problems tend 
to be, by their very nature, ill-conditioned, and any nostrum that patches the theory must surely be as likely 
to hurt as it is to improve the condition number; conversely, due to natural trade-offs that one might expect 
to have to make between theory and practice, anything done to improve the condition number must, it might 
seem, surely blemish the theory, in some sense or other. 

Surprisingly, upon considering several simple ideas and their implications, it turns out that the first more 
optimistic take on things is much closer to the final truth of the matter. One central point is that it is 
necessary, at least implicitly, to consider conceptual issues involving the linkage of the internal energy of the 
source density and the energy content of the exterior field. When this is done, it turns out that even the 
first-cut perspective mentioned above is, if anything, actually too pessimistic. For example, it turns out, 
upon closer examination, that all of the various relevant standard regularization techniques have a sound 
theoretical justification and, conversely, all of the considered theoretical nostrums lead to solutions which 
improve the condition number. When one steps back to consider what the implications of all of this might 
be for the application of the DIDACKS approach to other arenas, there is even more good news, in that 
the overall analysis leads to relevant strategies for handling various ill-conditioning issues that can arise in 
conjunction with Laplacian inverse source problems in these other areas. 

Conceptual Preview 

For any true measure of success over a broad range of relevant problem ares, a surprisingly large body of 
diverse ideas and corresponding techniques will have to be considered here. This led to inevitable difficulties 
in attempting to organize the underlying material. In particular, since there are various issues that are hard 
to discuss in a straightforward way using a completely linear sequential exposition, it seemed appropriate 
to set the stage in several different ways. Thus, many of the deeper issues subsequently encountered are 



raised in the remaining part of this "Preamble Section" and in the next "Overview Section," prior to the 
more standard "Introduction Section." Also difficult choices had to be made about what to leave out. In 
particular, since geoexploration and physical geodesy are the primary focus here and the material was chosen 
with an eye towards generality, some significant special topic issues are not addressed at all — for example, 
that of handling conductors in the field region for electrostatic inverse source problems or that of adapting 
the formalism to accommodate the magnetic dipole form (versus the electrostatic or gravitational form that 
is more easily dealt with in DIDACKS theory). In the end, the material that is contained in this article 
should be accessible to most of the readers from various other disciplines who might wish to try adapting 
it. It goes without saying that, in the end, the reader will have to determine what concepts are relevant and 
appropriate for his or her problem area of interest. In order to help orient general readers a brief survey of 
some of some of the main points made in the sequel will be given next. (The order followed here does not 
necessarily follow that of the discussion in the main body of the paper.) 

First, due to the fact that under the influence of gravity matter always attracts matter, the gravitational 
field energy is negative instead of positive. This is obvious since energy is released when separated bits of 
matter come together to form, say, some planetary or stellar body. (Conversely, the electrostatic energy 
density is positive and energy is required to assemble a charged body, since like charges repeal.) For less 
massive bodies, such as the Earth, internal stresses and/or pressures restrain the resulting static configuration 
so that it does not implode entirely in upon itself, but the final configuration may well be a quasi-static 
configuration, where energy transfer of some sort or other must be considered. (Obviously, tides, earthquakes 
and related terrestrial phenomena are manifestation of such energy transfer processes.) 

Next consider a concrete example of this gravitational body formation process. If a massive (but not 
too massive, so that black hole formation is avoided) isotropic isolated non-rotating cloud of gas coalesces 
into a spherical ball under the action of mutual attraction of its parts, then, at some point in time, internal 
pressure will restrain the configuration from further collapse. Even this simple example turns out to be 
fairly complicated, because time dependent thermal gradients and their effects on pressure must be taken 
into account. Thus the core of this configuration becomes heated as it initially compresses and this, in turn, 
influences the pressure of the core itself; moreover, subsequent radiative cooling will have a strong influence 
on the final density profile and on how soon a stable (or nearly stable) configuration is obtained. Detailed 
technical efforts to model planetary density configurations in terms of rotating gasses and incompressible 
fluids (as well as other forms of matter) have a long history and frame a pertinent part of astrophysics 
and geophysics [l6l [3]. For standard planetary bodies, these resulting axially symmetric configurations 
correspond to an ambient density profile that is homogenous at each depth (i.e., tangentially isotropic) and 
these resulting nominal configurations can be taken as defining a minimum energy or ground state density 
configuration. It is thus natural to remove this nominal reference field (which is labeled "normal gravity" by 
geophysicists [21) along with the included rotational effects and thus consider only density and gravitational 
field deviations from this equilibrium "ground state." Clearly these density differences from a nominal (or 
normal) configuration result are as likely to be negative as positive. The point being that these density 
differences indicate some sort of dislocation of matter that is associated with internal stresses and strains, 
so that they correspond to states of higher energy, which, in turn, means that the resulting external field 
differences also correspond to states of higher energy. It is thus clear, then, that when a suitable reference 
field is subtracted, and one is willing to express density estimates as excursions from nominal conditions, that 
the proper strategy corresponds to minimizing the external field energy, subject, of course, to matching the 
available data (in, say, a point- wise sense) . As noted in [12] this is precisely what DIDACKS point mass fits 
do: They are the collocation fits that minimize the overall (weighted) field energy subject to the constraint 
that the potential be matched at specified data points. [The DIDACKS approach also simultaneously 
minimizes the (weighted) field energy of the difference between some specified truth field and the field to be 
estimated, which, in itself gives a strong added motivation for subtracting off a nominal reference field.] 

As discussed in [12] , it is worth noting that this method of residual fitting, as it is referred to in [12] , is 



very useful within a general geophysical context. In particular, as pointed out at the end of the main body 
of that article (fT?), geophysical techniques generally always include subtracting off a standard ellipsoidal 
reference model from the total field and dealing with a quantity that is called gravity anomaly (or a related 
quantity called gravity disturbance) \2\. Often a further subtraction, which is called remove-and-restore (and 
that goes back to Fosberg, circa 1984), is carried out, which results in a localized field. At several other 
places in [12] the advantages of residual fitting techniques in DIDACKS based applications were also pointed 
out. 

Second, a degree variance analysis shows that when a spherical harmonic reference field of degree and order 
nine or higher is subtracted off, then the spherical DIDACKS weighted energy norm relationships for the part 
of the field that is left give almost identical results (i.e., norms) to unweighted energy norm expressions for 
spherical exteriors; consequently, when a suitable reference field is subtracted off, the interpretational issues 
that arise for spherical exteriors associated with the difference between weighted energy expressions and 
direct energy expressions can be ignored. It is also worth noting that this removal of a reference field has the 
effect of partitioning the density estimation problem into roughly two parts as well: (1) Core and deep mantle 
density estimation that is primarily tied to the chosen reference model chosen. (2) Geoexploration and other 
surface oriented geophysical density estimation areas that are primarily concerned with the remaining residual 
fields. The concepts and techniques presented here are probably of most direct interest in geoexploration 
and other near surface related problem areas, but they clearly can also be adapted to the core and deep 
mantle regimes as well. 

Third, while there is clearly much more analysis that can, and should, be done along similar lines — 
especially with regards to the connections of core and mantle density distributions, along with all of the other 
geophysical and geodynamical aspects — in Section[5] certain commonly used regularization criteria are shown 
to correspond to assuming a direct proportionality between gravitational field energy and matter dislocations 
and an analysis is presented which shows that this connection has a reasonable physical interpretation. Most 
of the other procedures that improve condition number correspond to assuming smoother density variations 
at the expense of choppy ones, which is one of the major overall themes of this paper. 

Finally, as an aside, although various types of point sources were considered in this and subsequent articles 
(including point dipoles and point quadrupoles), attention here is focused on theoretical and practical issues 
associated with mass density estimation from gravitational potential field information. Researchers in other 
disciplines that deal with inverse solutions to Laplace's (or Poison's) equation requiring dipole or other 
distributions as solutions can easily adapt the techniques presented here to their venue, so this is not nearly 
as restrictive as it may at first seem. 

Implementation Notes 

While, for most readers, it may be tempting to only ponder the various issues raised by this article, in 
order to have some true understanding of the issues raised, it is probably necessary to perform at least some 
numerical trial implementations. A basic DIDACKS point mass fit is very easy to perform. For example, 
in the half-space setting all that is necessary to do is to set up and solve d?]) using, say, the Householder 
triangulation routine in Lawson and Hanson [5]. Since this article discusses and tries to justify various 
implementation points primarily from a conceptual point of view, it may not be clear how simple and easy 
to implement many of the suggestion are. The goal of the remaining part of this section is thus to point out, 
by means of a few simple suggestions, how someone who is approaching these issues for the first time can 
gain some hands-on implementation experience in a relatively painless way. 

First, observe that since a DIDACKS point mass fit [i.e., V^^ specified by ([3])] satisfies the minimum 
collocation norm property, it is the potential with the smallest norm that matches the given (collocation 
point) data values of the specified potential W . This means, for example, that if Vjsi consists of a few shallow 
point masses only, then Vn —^ away from the specified collocation points, while if one adds more point 



masses and places them deeper the condition number will rapidly increase. A standard way of improving 
the condition number is to add a regularization factor that is a quadratic function of the parameter set to 
the cost function. For example, here the chosen modified cost function form might be 



<^=\\VM,-W\\l^+Tj2^k\ (1) 

fc=l 

where r is a small adjustable constant. Observe that this modified cost function only tends to exacerbate 
the tendency of Vn to "regress to the mean" since, for sizeable values of r it rapidly drives the values of nik 
to zero. Notice, however, if prior to performing the fit some given reference model has been subtracted off 
from W, then the tendency of Vn^ to "regress to the mean" implies that Vn^ will regress to values that are 
relative to some preselected reference model, which is generally a desirable property. Moreover, in the sequel 
it will be argued that under these circumstances one can conclude that this quadratic form is an appropriate 
expression for the potential energy that is associated with internal material dislocation. This subtraction of 
a reference model is called residual fitting in the sequel, just as it was in [10]. It should be a simple matter 
for the reader to test all of this out for him or herself. 

Suppose that no reference model is available for use in residual fitting: What then? Along simular lines 
that lead one to conclude that quadratic expressions like the one on the right hand side (RHS) of ([1]) are 
associated with minimum energy states, on can argue that, in general, the smoother the density profile, the 
lower the energy state. All of these considerations lead to a regularization factor of that is proportional to 

-^ ^ Yl ^iji'iTT'i ~ "^i)^j where the uJij are positive constants that are nonzero only for source points that 

are not separated too much. (Which clearly makes the overall density profile more or less smooth.) Again 
the interested reader can test this out directly using simple examples. 

Finally, it is worth noting explicitly that the DIDACKS formalism does not have a built in way of handling 
error measurements, like geophysical collocation does jGj, so some sort of preprocessing is necessary when 
measurement errors are present. As indicted in |11| . a certain amount of caution is generally necessary — 
especially when using raw data. To evaluate a candidate implementation, it is generally necessary to use 
some sort of realistic synthetic data and perform a Monte Carlo analysis to get a handle on the error behavior 
of the chosen implementation. (Even though standard covariance-based data preprocessing algorithms may 
give internal estimates of data statistics, they may not be completely reliable when, for example, downward 
continuation is present.) 



(ii) Overview 



Given the general historical and physical import of Laplace's and Poisson's equation, methods of solution 
for either one in some particular realm are of general interest since the methods employed may touch on 
solution techniques in many other problem areas. Thus while this article focuses on interior mass density 
reconstruction from given exterior gravity field information, many of the physical and mathematical strate- 
gies introduced here are quite general and can be either used directly (such as for gravity field modeling 
problems) or extended for use in other areas of applicability. The basic philosophy used here in approach- 
ing density estimation problems is energy minimization because physical systems in stable equilibrium are 
clearly minimum energy states, which means that an energy minimization based approach can be used to 
by-pass well-known theoretical issues of ill-posedness that are generally linked to Laplacian (or Poisson Ian) 
equation density estimation problems. 

This work builds directly on a mathematical framework presented in a previous article by the author 
that allows for the approximation of R harmonic fields in unbounded domains by (point) sources contained 



inside the (bounded) complimentary interior region [10] . At the heart of this approach is the idea of a 
Dirichlet-integral dual-access collocation-kernel (DIDACK) that has the ubiquitous form of the inverse dis- 
tance between some field variable point and a fixed source point. Because the field point and source point 
are assumed to be in disjoint regions this kernel form is bounded; furthermore, this restriction to disjoint 
domains circumvents theorems disallowing reproducing kernels that have this general form and, in fact, these 
kernels can be employed to obtain closed form expressions for energy norm inner products. Nevertheless, 
the resulting structure cannot be considered a reproducing kernel nor is it even a symmetric kernel since its 
two arguments do not share common domains. There are, however, natural connections of the associated 
space (DIDACKS) techniques, reproducing kernel Hilbert space (RKHS) techniques and especially geophys- 
ical collocation, which is a specialized R^ reproducing kernel technique [6]. The RKHS and GC connections 
to DIDACKS theory, as well as the general mathematical backdrop and various precedents were addressed 
in this previous article that also dealt with the R^ half-space, in addition to the spherical exterior setting 
[TU] . It also included an overview of relevant aspects of physical geodesy and a brief outline of the author's, 
as well as others, experience with point mass fitting. For the reader's convenience, the basic DIDACKS 
approach is reviewed in Section [2l but since no attempt is made in Section [2] to motivate or re-derive the 
basic DIDACKS mathematical relationships, in what follows it is assumed that the reader is familiar with 
the overall plan of approach. Connections of DIDACKS theory to various other mathematical approaches, 
such as the method of fundamental solutions, also exist and were addressed in a separate article in this 
series dealing with DIDACKS M^ and C theory [TT] and it is worth noting explicitly that many of the ideas 
developed here can be applied (or adapted) for use in these other mathematical settings. 

Underpinning the mathematical and physical basis of the approach to Laplacian inverse source theory 
presented here are the DIDACKS closed form expressions for gravitational field energy, which yield a con- 
sistent source estimation procedure and interpretation when supplemented with four realizations that are 
central to this presentation: 

1. The method of residual fits 

• subtracts off a nominal density profile so the results are as likely to positive as negative, 

• implies that any deviation from this nominal density profile are associated with an increase in 
field energy, 

• geographically localizes the source estimation problem so that the procedure can be readily 
adapted to the geometry and data at hand. 

2. When deviations from nominal conditions are under consideration, one can show that commonly used 
regularization procedures lead to a self-consistent physical interpretation and approach. 

3. Also conferring various implementation advantages, including greatly improved system condition num- 
bers, is the structured point source technique (SPST), where groupings of point masses are used (with 
each grouping often being some selection of nearby point masses, generally taken to lie on a regular 
grid) and all of the masses in each of the groupings have predetermined relative mass values so that 
only an overall mass scale factor for each of the groupings is determined. (Each of these groupings can 
also be regarded as defining a basis function.) This has the effect of replacing a given point source 
basis function with a distribution that is more uniform and spread out, not only in terms of its density 
representation, but also in terms of the effective potential and forces produced. As discussed below, 
one significant advantage of utilizing SPST basis functions is that a (family of) SPST basis function 
can be engineered to have characteristics that meet preselected requirements, in say the frequency or 
spatial domain. 

4. Finally, since the DIDACKS energy minimization approach is based on a cost function, it can be 
seamlessly integrated with other cost function based approaches due to the inherent additivity of all cost 



function based approaches. Moreover, alternative physical descriptions may easily be used to describe 
information content that was missing in the original cost function description: ^Tot = ^a + tb^b + ' • • , 
for tb > 0. {The point is that if the process described by $^(0?) is physically consistent with that 
described by $b(q?), for some (global) parameter set d?, then the minimum of <i>Tot(Q?) = minimum 
of $yi(a) = minimum of $B(a), while if the physical descriptions A and B are not consistent then 
T > should be chosen in such a way as to reflect the relative reliability of <i>A(Q?) and $b(q?) [in some 
instances the ratio of <J>a(q?) and $b(q?) (or its inverse) may be a direct expression of this relative 
reliability and suffices, as a general practical guide, for implementations]. 

The goal of this article is to articulate, clarify and amplify on these four points, while putting them within 
the overall context of a field energy minimization approach. Given that the exact nature of these various 
points, how they relate to each other and how they relate to field energy minimization may be unclear at 
this stage (especially with regards to items 2, 3 and 4), a few additional side comments are in order here. 

It has already been noted that items 1 and 2 are connected. Next, observe that items 2 and 4 are related 
since the general way regularization is added is by minimizing $(a) + rfl instead of $(a). Here r is an 
adjustable constant and fj is the so called "regularization function." (Usually an over bar is not used in 
denoting the regularization function fi, but an over bar is used here since the symbol O is reserved to represent 
the field region of interest) . Generally a regularization function (such as in Tikhinov regularization) is chosen 
solely for its condition number improving properties so t is chosen to be as small as possible, consistent with 
this overall goal. Here, however, the perspective is that ill-posedness is most likely a direct result of ignoring 
pertinent physical information about the underlying processes; consequently, the regularization process might 
be labeled "constitutive regularization" (versus Tikhinov regularization). (It is perhaps worth noting, that for 
some years alternative information based approaches to inverse source theory have been suggested [T4l[T5] V 
That is, as subsequently argued, the point is that energy is generally associated directly with internal source 
dislocations and that, from a general constitutive perspective, there is a direct correspondence between 
often used regularization forms and reasonable expressions for this constitutive energy; conversely, simple 
assumptions and a straightforward analysis of the nature of this constitutive energy leads to natural forms 
of regularization functions. It is also worth perhaps noting that a cost function based approach occasionally 
affords a easy means of collaboration. In particular, with regards to item 4, although experts from respective 
fields A and B may have knowledge of their own specialty only, they may be able to form a collaborative 
effort where each submits his or her own separate cost function for use in the final total cost function, yielding 
a unique composite optimal solution as a result. This point is germane since the formalism presented in 
Section 5, when implemented along the lines of the last example given in this section, should allow for the 
tight integration of seismology data and gravitational data. Finally, what is partitioned via 1 can also be 
added back by 4, so that there can be a subtle process of refinement of the total solution for certain relevant 
physical processes. 

Consideration of a few explicit examples may perhaps be necessary in order to clarify the main ideas 
behind item 3. Thus, first consider a source density region partitioned into a set of non-overlapping cubes. 
Although, in geoexploration, a collection of cubes (or parallelepipeds) such as this proves to be a very useful 
ansatz for gravity source estimation, there is an overall added level of implementation complexity due to 
the fact that the closed form potential (and gravity) expressions that result for each of the cubes [from 
integrating ([2|) below] are quite messy and this, in turn, clearly complicates any DIDACKS implementation 
(although such types of applications are not entirely unreasonable to consider — see the last part of Section 5, 
for example). When the RHS of ^ is taken over a cube (say) the resulting potential function on the 
left-hand side (LHS) can be regarded as simply defining "a cube basis function." On the other hand, one 
might consider a numerical approximation to this integral where the continuous distribution is replaced by 
a uniform grid of point masses, each of equal (but unknown) mass. A better way to think of this is simply 
to regard this grid of uniform masses, not as a numerical approximation to a continuous distribution, but as 
a distinct type of distribution, which is to say a structured point source (SPST) basis function (as indicated 
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above in item 3). In this case, a SPST basis function can be regarded as an approximation to a cube basis 
function and it can be translated from place to place, just like a cube can; although, obviously in this case, 
each SPST basis function must be indexed by a separate label indicating its location. Notice that this 
SPST basis function has only one undetermined parameter scale factor and that the total mass of the cube 
it approximates fixes this scale factor (or this overall scale factor can obviously be fixed by performing a 
DIDACKS fit — as discussed next). 

There are several additional points to be made here. First, DIDACKS linear equation sets for the source 
parameters that result from using SPST basis functions are exact, easy to implement and easy to solve. For 
example, when compared to a simple point mass fit, only the appropriate sums over the cubes have to be 
added to obtain the governing exact linear equations set. Second, the matrix size of the linear system to be 
solved is determined by the number of regions of the system to be modeled (i.e., number of cubes) and not 
the size of the basis function internal grid, because there is only one unknown source parameter per source 
region, or SPST basis function. Consequently, a very fine (internal) grid can be taken, if desired, without 
increasing the size of the linear system to be solved. Third, it is the dimensions of the cubes themselves (i.e., 
how close the source regions, or SPST basis functions, are together) and not the spacing of the internal grid 
that determines the condition number of the source parameter system. Thus the use of SPST basis functions 
corresponds to using a sort of "internal" or "structural" regularization and, as such, the SPST approach is 
directly related to other energy based regularization techniques (c.f., item 2). (The point here is that stable 
forms of solid or liquid matter have a certain inherent "energy of constitution" associated with either their 
molecular or crystal lattice structures and so long as this basic constitutive nature of matter is taken into 
account the actual amount of energy involved here does not matter since, aside from variations due to energy 
of deformation, the internal energies of material constitution are constant and hence ignorable so that only 
energy scales directly associated with deformations or dislocations need to be considered.) For example, 
instead of the collapsing cloud of gas or dust considered above, one might consider (a more realistic?) 
model of an assemblage of preassembled uniform density rocks or other objects, each of which, since it is a 
preassembled uniform clump of a given type of matter, has an innate internal energy of constitution. Under 
gravitation collapse, it is clear that, generally, the final state of such an assembly of matter will be in a 
minimum energy state, provided that there are insufficient pressures to cause excessive elastic deformations, 
and that under these circumstances the end configuration can generally be expected to be unique. Fourth, 
as hinted at already by this type of gravitation collapse example, the resulting software implementation of 
a SPST approach can be made very flexible so that it can be adapted to various shapes, sizes and types of 
objects (and corresponding source regions). 

This brings up the second type of SPST basis function implementation example, which shows that there 
is generally a connection between items 3 and 4. For the sake of concreteness, consider a case where there 
are three layers of unknown, but uniform density, and that each of these layers can be approximated by a 
single SPST basis function (each of which has an irregular boundary, in general). Further, assume that each 
of these slabs has an underlying common uniform point source grid, so that in matching the point sources 
to their corresponding basis functions it is only a matter of saying what grid point falls into what slab. It 
is then only a matter of deterring the overall (scaling) densities for each of the three slabs. Suppose that 
fl = {jii, fi2, A^a) represents the three source parameters (i.e., SPST basis scafing factors) of interest and that 
ff are a finite set of parameters that determine the location of the boundary surface between the first slab and 
the second slab and that all of the other surfaces are known (and thus fixed) by some other means. [The if 
might be, for example, representative surface points that determine a Junkins interpolation, or they might be 
surface spline points, or they might be, say, a set GC surface determination parameters (with the idea being 
that if the statistics of the surface height are known then these surface parameters can be used in (more or 
less) the same way that GC is used in performing geoid height estimation).] Then there are two sub points: 
first, one can simply perform a DIDACKS procedure to determine the total parameter set (a.) = {jl \ ff)'^ 
by minimizing $(a) via standard nonlinear least-squares (NLLSQ) optimization means. Here, for example 



if variations in the surface are due to mechanical stresses in the slab (associated, for example, with flexure 
of the slab itself) then a direct energy cost function can be associated with these surface parameters and, in 
accord with item 4 (as well as item 2), an energy regularization cost function can be added to $(d), thereby 
both improving the estimate and its numerical tractability. Second, in accord with item 4, additional data 
source types for determination of (a) can be entertained besides gravity data, provided they can be written 
directly as a cost form process (here one obvious example for geoexploration problems might be to include 
seismology data). 



1 Introduction 

First, the overall mathematical setting can be described in terms of the Earth's gravitational field Ge{X) 
over the unbounded exterior domain, fJ C M' , where the mass of the Earth is assumed to be contained 
inside the compliment of fl, denoted fl' . This vector field is derivable from a scalar field or potential We{X): 
Ge{X) — WWe , where a positive sign on the RHS, rather than a negative one, is used to conform to the usual 
convention adopted in physical geodesy [2] (other conventions were noted in [10] , but are of no immediate 
concern here). Given this assumed default linkage of vector field and scalar field, only scalar potentials and 
their sources need be considered in the sequel. The Earth's gravity field arises from some source density 
Pe{X') that is contained inside a bounded source region, fl'g C fJ' C K . (Primed variables will generally be 
contained in the source region and unprimed ones in the exterior region, so a prime has been affixed to the 
source region symbol.) Thus Poisson's equation, V^ We — —4:npE, holds for the whole of R and Laplace's 
equation holds for SI since Pe := there: V^ WeIX) = for X £ fl d M'^. The potential field in Q and its 
density in Q'g are then linked by: 

which is the integral form of Poisson's equation. [The question of particular units to be chosen is by-passed 
here, so a constant factor may need to be inserted on the RHS of ^.] 

Within this overall mathematical context there are two significant broad historical areas of research to 
contend with: (A) The determination of the Earth's global density profile, which along with the determination 
of the Earth's shape can be said to comprise the central issues of physical geodesy [2]. (B) Problems associated 
with density determination for more localized regions arising form petroleum and mineral geoexploration 
efforts. Here (A) addresses either deep mass distributions or shallower densities that do not vary abruptly, 
while (B) deals with near surface densities, and regions of abrupt change are often of special interest. In 
some sense, geodynamics [IZIIT] addresses issues that span both of these scales since it deals with phenomena 
such as plate tectonics and earthquakes, but since only configurations in static equilibrium (i.e., non-time 
dependent ones) will be explicitly considered, these problem arenas are not be addressed here in any detail. 
Note that Pe and We can be conveniently partitioned into two parts corresponding to (A) and (B): Pe = 
Pa + Pb and We = W^ + Wg and that this partition simplifies a host of related interpretational issues. 

Here (A) goes back to the origins of potential theory itself and already had a long associated history by 
1900 [in]. The realization that the problem of attempting to estimate interior mass density profiles from 
exterior gravitational fields is ill-posed goes back to Newton himself who showed that a uniform spherical 
mass shell and a point mass at the center of this shell produce the same exterior field provided they both 
have the same total mass. Here it is assumed that the relevant aspects of this part (A) global field can 
be captured in terms of spherical harmonic expansions. Recent progress in this area has been spurred by 
deployment of advanced satellite systems, such as the ongoing dual satellite Gravity Recovery and Climate 



Experiment (GRACE) (as well as the GROCE mission) [2]. As previously noted, when spherical harmonic 
expansions of degree and order 9 (or higher) are subtracted off, then the fitting results obtained from the 
spherical weighted energy norm (i.e., the integral norm) match those of the spherical energy norm very 
closely. Furthermore, if some care is exercised, then when a spherical harmonic field of degree and order 120 
(or higher) is used as a reference and removed then the half-space energy norm can be used. 

2 Synopsis of DIDACKS Approach 

This section briefly outlines the basic mathematical formalism developed in [10., where the main focus was 
point source Laplacian field reconstruction problems (i.e., gravity field modeling and estimation problems). 
First, consider the general DIDACKS plan of approach. This approach is based on minimizing energy 
based norms of the difference between some point mass (or more general point source) model potential v{X) 
and some given canonical reference (or truth) potential wlX). This can be restated directly in terms of 
minimizing some cost function $' = ||w — 1(;||'^, where || • || is the (possibly weighted) energy norm of interest. 
For point mass basis functions the potential model becomes 

^(x) = Gy^;^^. (3) 

In ([3]) G is the Newtonian gravitational constant « 6.672 x 10~^^m'^s^^Kg^^ 2 . Here it useful to introduce 
scaled versions of the potential functions and to denote them by capitol letters so that V = v/G and 
W — w/G so that the factor of G need not be considered in the sequel. The relevant cost function thus 
becomes (where $ :— <I>'/G): 

^^\\V -W\\^ ^\\vf -2{V,W) + \\W\t (4) 

In the DIDACKS approach, since (weighted) energy norms for field energy expressions are used, the problem 
becomes to minimize $ — HV' — W^|||;. := HJn Mjl^^~^W^P (PX, where /ij = ^J'{X) is the weight function 
(which may be set to one). Besides this norm, (weighted) energy inner products will also be needed in the 
sequel: (F, W)e, := HL, ^'o^'^ ' ^W(fX. 

General Mathematical Considerations 

The notation conventions of [1^ are followed here. Cartesian coordinates are used in the sequel: X = 
{x, y, z)'^ G R^ and arrows are used to denote R'^ vectors, while n dimensional vectors (for n > 3) are denoted 
by lower case bold letters and their associated multidimensional matrices are denoted by upper case bold 
letters. Further, Rq will denote the radius of the sphere associated with VLq and the coordinate origin will 
be chosen to coincide with the center of this sphere so that fio = {X G M'^ | \X\ > Rq}- Likewise for the 
half-space case, the origin will be chosen to be in the plane dQ.i and the positive z-axis will be taken normal 
to the plane so that fii = {X € M'^ | z > 0}. These two geometries will be denoted Q.j, where j = or 1. 
(Observe here that the general visual shape of the subscript matches the shape of dVlj itself.) One other 
aspect of DIDACKS theory is worth noting, before addressing mathematical preliminaries. 

Here general relationships that hold for both geometries of interest (fio and Vti) will be considered. For 
X e Vtj (j = or 1) consider a vector field, G{X), derivable from a scalar field W{X): G{X) = WW, where 
X g K^ and all the sources are assumed to lie in some bounded "source" region, fl'g, C J7'-. We restrict 
ourselves from now on to potential functions that fall off at least as fast as 1/r as r ^ cx3 in f2j. 
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R^ Half-space (Qi) 

For concreteness, first consider the minimization process in M lialf-space, f2i. Here z > characterizes 
the region of interest (rJi). It is clear that we wish to minimize a quantity with the following energy- like 
form (where the factor of Stt has been inserted since it occurs in the electrostatic field energy expression): 

\\V -W\\l^ := 1- JJJ \VV -VW\' d'X . (5) 

2 2 2 

Here in general \\V — W\\„ := II^IL + II^IL ^ 2(y, W)„ , where, of course, the energy inner-product 
expression introduced after (j4]) is to be used for (V, W)„ . 

In particular ii V ~ V^^ is a point mass model of interest (with N^ masses) and W is an appropriate 
given field, then Vn^ = J2k=i mk/\X - Xj^\ where X e ^i and X'^ e Q'g^. Further if Ik := \X - X'^\, then 

\\Vn.-W\\1^ := \\W\\l-2j2mk{i^\ T^)^^ + ^ ^ m, mfe,(4-i,£-,i)^^ . (6) 

Observe that the first term on the RHS of ([B|) is a constant term. Taking the partial of Equation © with 
respect to mk" for k" = 1, 2, 3, . . . , Nk and dividing by two yields a linear equation set that can be easily 
inverted for the m^ values, provided the required inner products can be easily computed: 

^ Tk,k'mk' = Ak, (7) 

fe'=i 

where n,k' = {ik\^k'\ and A, = {WJ-\ . 

The DIDACKS formalism allows for the explicit closed-form evaluations of all the inner products occuring 
in ([T])- In particular, the energy inner product in this case is 

and, as shown in [TU], 

ii-\W)^^=W{x'k,yl~z'k)/4, (9) 

which can be used to evaluate the inner product terms Tk.k' and Ak in ([7]): 

M.^ Spherical Exterior (fio) 

Here flo ~ {X £ K^ | \X\ > Rq} describes the region of interest; however, matters are more complicated 
than they were for the half-space case. First consider two general admissible functions / and g (that is 
functions that are harmonic in f2o and which tail off to infinity at least as fast as 1/r). The energy inner 
product in this case is 

(f^9)^.---iffhf-^9d'X. (11) 
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{f^g)^:^-^[U,{frg)da^-^[[[Dr{frg)\ da, where D. := |- (12) 



The inner product for the "integral norm" [10] is also very useful here: 

4:TT J J ' ^'' ^' ' 4:TT J J '■ ' '"' ""^ r=Ro " ' - - ' Q^ 

and the RHS of (fT2| follows from the evaluation convention given by 

11 f{r,9,4,)da:^ J J [f{r, 9, cf>)]l^^^sm{9)d0 d^ (13) 

for standard polar coordinates r, 9, (p. Here, as in [2], a and da have the standard meaning when associated 
with the integrand J{X). Likewise, let the surface inner product be defined as 



(«, w)^ := (l/47r) jj v{r,9, (t))w{r,9, cf)da. 

With these definitions it is fairly easy to show [10] that 

(/,.g)i = 4i?o(/,5)E„-^o(/-5), ■ (14) 

For DIDACKS applications over R"^ spherical exteriors the integral norm is more important than the 
energy norm since closed form inner products can easily be computed from the following expression 



(/, g\ ^ PkW [Pk] , (15) 

where Pk ~ \Pk\, with 



R? 



^fc = ^TTTT U^ (16) 



k 

k 
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for some point mass location X'^, . Here, the integral norm can be reinterpreted as a weighted energy expression 

m 

{f,g\^^jjjr-'Vf-^gd^X (17) 

no 
so that jiQ = Ro/r {r :— \X\) is the associated weighting factor. 

3 Inverse Source Theory Prologue 

Petroleum and mineral and geoexploration are ongoing and historically significant research areas, where a 
considerable amount of time and effort has gone into exploring various alternative approaches and there 
is an extensive associated literature. When a source distribution of interest produces a well delineated 
signal that can easily be separated from the background distribution, it is possible to simply compare the 
resulting potential pattern with some precomputed one. Historically, this "forward solution" technique has 
been (and is) popular and it was probably the one first used [T]. However, in the literature, when more 
sophisticated approaches are called for, the issue of the proper gravitational source estimation algorithm 
to use immediately becomes less clear. While most of the approaches seem to work, all the currently used 
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sophisticated approaches in this area entail a certain amount of arbitrariness or lack of physical motivation, 
which seems to be inherent in the foundations of all the approaches. 

Conceptually, in terms of the V and W above, one might frame the ideal goal as being to directly minimize 

$ = JJJ^, {pv ~ PwY d?X , where pv and pw are the modeled and reference source terms. However, when 

^ . . . 

only W and its derivatives are known in Vlj , there is no apparent way to effectively frame this minimization 

goal in a workable fashion. This should be obvious from the fact that a continuous density like pw is generally 
nonunique, which, in turn, is clear from the fact that all spheres in a given region that have the same total 
mass and that share the same center produce the same external field — regardless of their radius. As noted 
above, what is generally overlooked, however, is the fact if one takes into account the internal energy of the 
medium that is associated with stresses and other physical processes a local minimization of total energy will 
(generally) result in a physically unique situation (i.e., density estimate), because all physical problems in a 
static-stable equilibrium have a energy minimum underpinning. Those situations that do not have a unique 
energy minimum are of special geophysical interest since they generally represent earthquakes, tides, core 
rotation or other geophysical situations, where dynamics (and the energy forms associated with it) must be 
considered. All of this was discussed at some length in the first two sections of this article [Section (i) and 
Section (ii)]. The solution to the foundations of inverse source problems proposed there involves utilizing 
the DIDACKS approach (to account for external field energy differences) in conjunction with augmented 
energy-like information added as a regularizing factor. Residual fitting also pays a central role in the physical 
interpretation. Energy as a basis for studying earthquakes has been proposed by others. The material stress 
energy models in these studies is often very detailed and goes far beyond the scope of what can be included 
here, but integrating the DIDACKS approach with these considerations is clearly an avenue that warrants 
future effort since the contribution of external gravity field energy has been generally ignored in this arena. 
In the current paper a constitutive regularization approach is taken and the goals are much more technically 
modest. The goal is simply to physically justify easily implemented internal energy minimization approaches, 
where flexibility and ease of implementation are maintained as a primary goals. 

It is worth explicitly noting that the general role of energy minimization obviously has not gone unnoticed 
historically. In particular, Kellogg [4J explicitly points out the role or field energy minimization in electro- 
statics via Dirichlet's integral [H p. 279], but there are clearly many other historical connections that can 
be pointed out in this context. It is also worth noting that when the foundations of geophysical collocation 
were debated by Krarup and Moritz, Moritz put forth a statistical interpretation (which is commonly called 
least squares collocation) that eventually won out, but that Krarup put fort the idea of a weighted energy 
minimization approach based, effectively, on the RHS of (fT7|) . As discussed in the last section of [10], since 
the goal was to give an interpretation to GC his idea was to use the Krarup kernel (rather than l/£fc used by 
DIDACKS theory) and to argue for the physical importance of energy minimization. Part of this debate can 
be glimpsed from some side comments in early geophysical colloquium proceedings. As pointed out in }10j . 
DIDACKS theory turns things around conceptually and abandons the pretext of a symmetric reproducing 
kernel (SRK), which largely severs the direct connections to geophysical collocation, while it keeps energy 
minimization and the fundamental of the form \X — X'^\~^ in tact. In the end, as argued in this paper, this 
also has the effect of keeping connections open between density estimation and energy minimization. 

4 DIDACKS Implementation Discussion 

As pointed out in Section [51 there are two field regions of interest, the exterior of a sphere and positive 
half-space, denoted by fij (for j = Q and 1), respectively. The corresponding energy norm for these two 
regions can then be simultaneously referred to as || • || . The corresponding DIDACKS norms for these two 
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regions can likewise be referred to as || ■ |L where 

MId, -II ■ IIei a^d II • llj^^ :-/c|| • llj . 

Here A: is a constant which can be chosen to preserve connections to units of energy [i.e., k = l/(4i?o)]; 
however, since the resulting values of source parameters for a DIDACKS fit do not depend on the specific 
value of this constant, it is generally more convenient to simply set k = 1. For fli, \\ ■ |L and || • |L are the 
same so the question of which to use does not arise; however, for f2o it would appear that for inverse source 
problems there is some question as to whether || • || or || • || should be employed as the major tool of 

analysis; where, of course, || • ||g has a direct bearing on assorted energy based arguments, but || • ||j-j is more 
mathematically amenable. Nevertheless, as previously pointed out, this dichotomy goes away for geophysical 
problems if a suitable low degree and order spherical harmonic reference is subtracted off (and then restored 
at the appropriate time), because the resulting residual field has no low degree and order content and the 
two norms, in this case, are nearly proportional. That this is so can be seen from a full spectral analysis. 
That this is so can also be seen by taking stock of p7|) . where it was already observed that fiQ — Ro/r is 
the effective weighting factor for the weighted energy norm in this case. The point here being that after a 
reference field is subtracted off, the remaining residual field attenuates very rapidly as r increases so that 
only values of the field close to Rq make significant contributions to the norm and in this "near field" region 
fiQ is approximately constant. 

There are two direct consequences of the DIDACKS (weighted) energy minimization approach: 

2 

(A) Since \\V ~ W\\„ is minimized, the resulting fit will be the one which minimizes the (weighted) energy 
difference of the error field [which by definition has the potential form (V — W)]. 

2 

(B) Since the GC property is satisfied the resulting field will be the one which also minimizes || V|L subject 
to the constraint that the sample field data points be matched (which is the replication property, so 
that, for example, for point mass fits V{Pk) = W{Pk) at all the specified data points {Pk}kJ^i )• 

While the interpretation of all this is all relatively straightforward, there are various issues that warrant 
consideration and further clarification. The first issue to be reconsidered is the sign of the gravitational 
energy itself. 

Negative Gravitational Field Energy and Residual Fitting 

Amplifying slightly on the discussions in previous sections, the electrostatic case will be compared with 
the gravitation case. Both electrostatics and gravitation are inverse square law forces and in both cases, for 
source free regions, the associated scalar potentials obey Laplace's equation. In the electrostatic case, the 
forces between two charged bodies are proportional to the product of their charges, while in the gravitational 
case the forces are proportional to the product of their two masses. Aside from the fact that masses are 
always positive and charges are not, which causes some minor interpretational issues here, there is one 
fact that cannot be ignored: like charges respell and like masses attract, so gravitational forces are always 
attractive. (The sign differences associated with the choice of the gradient of the potential may complicate the 
identification of gravitational potential and potential energy, but this is not a issue that needs to be addressed 
here.) The fact that gravity is always attractive means that the gravitational field energy is inherently 
negative, unlike the electrostatic case. To see this consider what happens when a set of gravitational or 
electrostatic sources are assembled from an initial configuration that is well separated (i.e., out at infinity), 
which is, in general, how one computes the field energy. In the electrostatic case, when a collection of like 
sources are assembled it is clear that positive work must be done to overcome the mutual repulsion of the 
charges. From well-known arguments found in standard physics texts the resulting electrostatic field energy 
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has the form = constant x -g^ JJ E ■ E cPX > 0. However in the gravitational case work is released during 
the assembly process (which, as noted in Section (i), is invariably gravitational collapse process) so that the 
work of assembly is thus — g^ // G • G cPX < 0. 

2 

Clearly the fact that the gravitational field energy is negative and thus related to the negative of ||W^|L 
is of prime interest here, since this difference in sign, at a minimum, is somewhat unsettling from the require- 
ments of a consistent physical interpretation. First, however, notice that, as discussed in [10', conditions 
(A) and (B) above, by themselves, provide sufficient direct motivation for handling modeling and estimation 
problems, since these conditions imply that the gravity error field is a minimum [by condition (A)] and that 
the fit is the most conservative one consistent with the given (point) data set [by condition (B)]. With regards 
to both conditions (A) and (B), notice that since the absolute value of the field energy is minimized, smaller 
(and therefore more conservative) overall excursions are emphasized over larger excursions. Also observe 
that (A) implies that it is not the absolute value of the field energy itself that is minimized, but the positive 
absolute value of the energy of the error field itself and this is clearly desirable. In general, as previously 
noted, residual fitting is a part of the DIDACKS approach, which means that W usually represents not the 
raw gravity field itself, but a residual field where some suitable well defined base reference function has been 
subtracted off. The resulting residual field can then be assumed to be zero-mean in the sense that when a 
point mass fit is performed on it the resulting point mass values occur in roughly equal positive and negative 
proportions, which is to say that X)fc=i "^fe ~ 0. As pointed out in Sections (i) and (ii), and backed up by 
analysis in Section \5\ this means that positive field energy is associated with (nonzero) source excursions, 
which provides a direct explanation of the negative gravitational field energy mystery noted above. This also 
clearly provides a direct explanation for the possibility of negative mass values arising in conjunction with 
point mass fits. Condition (B) thus must really be considered as holding for a residual field, where all the 
field excursions are to be considered as excursions from zero, so that the energy of the difference fields can 
always be considered positive. This means that the absolute value of the residual field energy is minimized 
subject replication constraints, which is clearly desirable. 

The Point Source Support Problem (PSSP) 

Clearly for standard point mass fits used in gravity modeling or estimation the DIDACKS approach does 
not generally require numerical discretization or numerical integration since the underlying linear equations 
are in a closed form; however, in this case ease of software implementation does not necessarily translate 
into uniformly care- free applicability. For a straightforward point mass fit, depth and spacing issues must be 
handled more-or-less correctly in order to obtain acceptable results. One common mistake made in utilizing 
point masses for modeling purposes is to not place them deep enough. Consider a point mass model based 
on linear equation ([7]). From the form of Tk^k' specified by (jlOp [where the overset tilde indicates the use 
of normalized basis functions] it is clear that as point masses are moved closer to the surface or further 
apart horizontally they clearly become less correlated (that is, Tk^k' becomes smaller). In the limit that the 
point masses are all near the surface they match the prescribed potential values at the specified locations, 
but the given potential model itself falls off to zero very quickly at locations away from those specified 
data points. In fact, under these near-surface circumstances a single point mass fit to a single data point 
behaves very much like a Dirac delta function. This is clearly an undesirable situation and to overcome it 
one must place the point masses at a fairly sizeable depth. Alternatively, as the point masses are moved 
closer together or placed deeper they quickly become overly correlated {Tk^k' — > 1), which, in itself, can lead 
to wild and unexpected variations in the resulting field model at points away from the prescribed reference 
field data points. These results clearly hold for the region fig as well and they hold for other point source 
types, such as dipoles and quadrupoles, as well. Often a fairly fine line between these two just outlined 
opposing and unwanted behaviors must be negotiated. A good balance of spacing and depth must be struck 
and when potential data locations of various heights and spacing is involved, it can become a very difficult 
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(or even nearly impossible) problem to overcome. This problem is thus labeled the Point Source Support 
Problem (PSSP), because it has to do with support issues associated with the underlying point source basis 
functions themselves. Clearly these problems are dependent on the choice of basis function, which is one 
reason that SPST basis functions were introduced in Section (ii). As previously discussed, various other 
techniques can also be utilized to overcome these issues and much of the rest of this article will explain 
relevant aspects of them within one context or another. The main means for overcoming the PSSP discussed 
in this section are residual fitting and spectral bandwidth decomposition. Various types of regularization 
and various distributed types of sources will be considered in the context of inverse source determination 
theory. Condition number considerations are closely tied to the PSSP and are also discussed below. 

Residual Fitting and Spectral bandwidth Techniques 

Next consider residual fitting as applied to gravity estimation, gravity modeling or source estimation, 
where it should be normally considered as an integral part of the DIDACKS approach in these areas. There 
are three primary reasons that the residual fitting technique is so effective for point source problems. First a 
certain number of degrees of freedom are always tied up in reproducing the general trends of the underlying 
reference model and when these reference trends are no longer present these additional degrees of freedom are 
freed up and can be used for modeling finer detail. Since simultaneously fitting a rapidly changing gravity 
field (which tends to require shallow point sources) and a field with long term trend properties (which tends 
to require deep source placement) is often difficult at best, residual fitting can be used to eliminate much (or 
most) of the long term attributes to be fit so that the regional part of the fit becomes, not only much more 
accurate, but easier to effectively implement. (Thus helping to overcome PSSP issues.) The second reason 
residual fitting is effective is associated with the collocation replication property, which DIDACKS fits obey. 
As previously noted, for techniques satisfying the collocation property, the fit usually digresses to zero away 
from the specified (field) data points; however, when a reference is subtracted off, this natural digression 
will be to the underlying reference model itself so that there is a natural attenuation built in. (This, in 
itself, clearly also helps alleviate point source support problems.) The third reason will be addressed next 
by itself and has to do with consistency of physical interpretation and is tied to more general field energy 
considerations, as previously discussed. 

After a residual fit has been performed to model the field to a certain physical scale, the entire process 
can be repeated and when such a series of residual fits is performed there is a synergistic effect. First since 
the (weighted) energy norms tend to fit the longest wave lengths first, the first fits performed (naturally, 
with sources chosen to be at a greater depth) will account for that part of the field that tails off more slowly 
with altitude. In turn, when this part of the field is treated as a reference and removed only the shorter wave 
length and more regionalized part of the field remains to be fit. The whole process can then be repeated 
as needed. In conjunction which this repeated residual fitting process note that it is important to remove 
the longer wave lengths present at each of the successive stages, or much of the error at each stage will be 
folded into the parts of the field that are be to fit subsequently. In this connection, it is worth explicitly 
noting that for the integral norm a degree variance analysis (or harmonic Fourier series analysis for the 
half-space energy norm) shows explicitly that a strong premium is placed on correctly matching any longer 
wavelengths that happen to be present. There is one negative aspect of residual fitting techniques. Since 
residual fitting techniques work primarily due to preconditioning of the 'signal' (so that it can accurately 
be fit by point sources), in general a good resulting point source fit will entail sources that are deeper than 
one might normally expect. This, in turn, leads to associated condition numbers that are large. If signal 
errors are present or if source estimation is the main goal, rather than modeling, clearly there may reason 
for concern. 

In many cases a spectral bandwidth approach can be used in place of residual fitting. With regards to 
the global part of the field a spectral bandwidth approach simply entails dividing up a spherical harmonic 
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expansion of W into various degree ranges so that each resulting spectral band has well defined physical 
attributes. Since such spherical harmonic reference fields are both accurate and readily available, it is 
assumed, that at the very minimum, that some low degree and order reference field will be subtracted off 
and used as a base reference for either a residual fit or spectral band approach. For local or regional fits 
where data is specified at some survey altitude, one might use digital filtering techniques or a fast Fourier 
transform in order to obtain various spectral or frequency bands. With sufficient ingenuity, interested readers 
should be able to figure out any required further details for such implementations so they are not discussed 
here. There are, however, two further points that are worth commenting on in this subsection. 

Considering the spherical case, it is natural to assume that the surface of the earth is to be taken as 
coincident with flj; however, when the data is specified at a fixed altitude pB]) fixes the point source depth 
at corresponding depths that may be totally inappropriate for the associated frequency content. Clearly one 
solution is to simply upward (or downward) continue the original data by using GC to an altitude that will 
produce an acceptable point source depth for the required spacing. In practice an alternative technique that 
generally works quit well is more appropriate. For concreteness, consider a DIDACKS point source fitting 
problem based on the spherical exterior geometry using the spectral-bandwidth approach. Further, suppose 
that data is given on the surface of the Earth, which is specified by \X\ — Re- (In practice, for a regional fit it 
is natural to take the origin to be directly under the center of the region along the ellipsoidal normal direction 
and at a distance that best captures the ellipsoidal curvature effects over the total region of interest.) Next 
determine what the appropriate spacing and depth should be for an ideal fit. The approach is then simply 
to consider Rq to be a variable (< Re) that is to be set to a value that will insure that P^ produces this 
desired depth — for data sampled at the correct spacing. This technique entails no loss of consistency since 
these deeper point sources are associated with a field region that may, in fact, be naturally taken to have a 
much smaller Rq. In particular, there is no reason to ignore the (weighted) field energy between Rq and Re 
produced by these deeper sources by arbitrarily taking Rq — Re- (In any case, for deep sources the source 
exterior region clearly has a boundary that is somewhat below Re and there is no real reason for thinking 
that Rq = Re is the correct boundary for weighted exterior field energy minimization over this part of the 
field.) This same technique can also easily adapted to the geometry specified by Qi. Here it is a simple 
matter to move the plane dili deeper, which forces the associated point sources themselves deeper. 

Finally, in the present context, it is interesting to note that a simular boundary adjustment technique 
to that just described can be used to produce a norm that minimizes energy over a region bounded by 
two planes (or weighted energy over a region between two concentric spheres). Although the technique is 
general it is easiest to describe it in terms of a single point source for the geometry Hq. Thus, suppose that 
fig = {X e M^ I \X\ > Ro} and OJ = {X e K^ | \X\ > Rq} describe two spherical DIDACKS regions with 
the same origin and Rq > Rq. Let Sl@ := flQ — Hq be the region of interest, then from (fT5|l there results: 

T)[w, q\ ^iQ, n^] = \Pk\ w{Pk) - \P*,\ w{Pt) (is) 

with Pk from (1161) and 

(19) 

Clearly a like expression follows for the region bounded by two parallel planes. 

Condition Number Considerations 

Commonly available singular value decomposition (SVD) or Householder triangulation routines are ap- 
propriate for solving the DIDACKS point source determination linear equation sets. (Generally, the amount 
of processor execution time is so minimal that it is not a real consideration and, thus, except for rare cir- 
cumstances, the universal reliability of the solution from a SVD or Householder triangulation algorithm is of 
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much greater importance than execution time.) The data repHcation property allows one to verify that the 
implementation has been performed correctly. The eigenvalues, which should be all positive, can be readily 
obtained from SVD routines. (Here it is worth noting that for those few cases where one might require 
ruk > linear inequality constraint software can be employed, but note that one should generally validate 
the output of this software [9 .) The condition number (C#) is taken here to be the ratio of the largest to 
the smallest eigenvalue of the T matrix in this linear equation set. 

As discussed above, for modeling problems, when the masses are too shallow the C^^^ will be too small. 
When the masses are too deep the C^^^ will be too large. For source determination problems matters are 
somewhat different and the C^ should be somewhat smaller. The point is that a large condition number 
is generally associated with large variations (and thus uncertainty) in the estimated masses and this is 
obviously associated with uncertainty in the prediction of the mass density itself; moreover, a large C^ 
indicates that any measurement errors will tend to be magnified by a like amount. Special techniques that 
lower the condition number, while simultaneously overcoming the point mass support problem, have been 
emphasized elsewhere in this article (obviously, many of these techniques are also appropriate for use in 
gravity estimation and gravity modeling problems as well). A guiding principle is that the less certain one 
is about the fitting results the lower the C^ should be; moreover, for modeling problems in general one does 
not care about the mass values themselves — only the results. This means that unphysical mass values are 
perfectly acceptable if they produce a good fit. In this context, it is perhaps worth noting that for low degree 
and order global point mass NLLSQ fits good results are associated with C^ ^ 10^*^, but for point mass 
modeling with various fixed locations, one would generally expect somewhat smaller condition numbers than 
this threshold (results for grids over interior regions were noted in |12| . but the same gridding techniques 
can obviously be profitably used for exterior regions as well). 

Before proceeding with the discussion of other fitting techniques in this next section, a word of caution 
is in order. Some sort of experience with point mass fitting is probably required before attacking real world 
estimation or inverse source determination problems. Thus, it is suggested that the reader interested in 
these areas gain as much experience through synthetic modeling as possible by working with various field 
models W, which are chosen to have properties that are as realistic as possible. For estimation problems, 
such modeling allows one to check the produced field values by the intended spacing of the point set at 
various locations away from the field sampling points. Likewise, for source estimation problems one can test 
the predicted source values against the ones used to produce W. 

Finally, it is also worth noting that the PSSP can be overcome generally by inputting a tolerance to 
the linear inversion software so that unwanted small eigenvalues are ignored (this is generally a very strong 
form of regularization) , thus allowing a very tight point source grid spacing while preventing large source 
values. Here is also worth noting that when normalized basis functions are used, the largest eigenvalue for 
the system (O, or its spherical analog, is obviously bounded from above by Nk- 

5 Inverse Source Determination Techniques 

For source determination problems C:/^ concerns must be addressed and there are two primary means of doing 
this: regularization techniques and basis function modification techniques. As just discussed in Section [H 
for optimal fits, residual fitting techniques generally overcome the PSSP, but generally at the expense of 
large C#'s; however, residual fitting is still an important source determination technique, since it works syn- 
ergistically with basis function modification techniques and regularization. The primary types of alternative 
basis functions to point sources that will be considered are structured point sources (which consist of an 
aggregate of point sources). Basis functions that yield continuous source estimates are also considered. In 
particular, since continuous sources are inherently nonunique, the primary tool considered in this regard is 
parameterized continuous source estimation. Other parameterized continuous source techniques have long 
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been used in geophysical inverse source theory [5] 

Regularization Techniques 

Before proceeding to the analysis of the physical basis of specific regularization approaches, it is useful 
to recap part of the analysis given in previous sections from a slightly different perspective. 

As noted in Section [3l for density estimation problems one would ideally like to minimize an expression 
like $ = JJJc,! {p^ — P^Y d^-^ where p^ and p-^ are the modeled and reference source terms and fig. is 

the source region associated with the two field regions of interest (fij); however, there appears to be no way 
to mathematically frame this minimization goal when only W and its derivatives are known in fij. It is 
easy enough to see that this must be the case from purely information content alone since W is harmonic 
in f7j and thus is specified by its values on d^j, whereas p-^ has many more degrees of freedom and is 
not determined by its values on a surface. As noted in [10], for a chosen closed surface in J7j, DIDACKS 
theory links specification of potential values to source values specified on a corresponding closed surface in 
the source region. From this observation and the fact that equal mass concentric shells produce the same 
external field, one can conclude that much of the ill-posedness of the density problem are associated with 
source depth issues. Two aspects, however, are clear: (1) Regularization, in any reasonable form, should 
help to stabilize the source estimates and thus generally provide more reliable estimates. (2) Gravitational 
inverse source problems minimize energy in some sense or another since all physical systems in static-stable 
equilibrium have energy minimum underpinnings. The general thrust here is thus to try to link these two 
aspects in the DIDACKS approach to source determination. As previously noted, residual fitting is linked 
to the interplay of these two aspects. 

In order to ascertain some of the underlying issues involved, consider a straightforward application of the 
point mass fitting theory presented in previous sections (without regularization). In this case the approach 

2 2 

is based on the minimization of || V^— iy|| ~ ||^— W^Hp, for held information specified in the region ilj (for 
j = and 1). (In what follows the norm expressions will be written in terms of energy norms.) As before 
it is assumed that an appropriate reference function has been removed from the specified function W prior 
to fitting, which entails the removal of a reference density from p^ as well (but this removal may only be 
implicit). Residual fitting helps here since residual fields have reduced low frequency content, which allows 
for shallower point masses placements of the remaining sources. The point being that the ill-conditioning 
arising from the source depth ill-posedness mentioned above can be overcome by introducing a natural source 
depth stratification. This helps to control one cause of innate ill-conditioning, but there is another one that is 
associated with how close the sources are together. Thus, in general, the condition number will still be much 
too high as the grid spacing becomes small. As previously noted, a large condition number is unacceptable 
in this case for three reasons: (1) The mass values will tend to vary wildly from one point mass location to 
another, so a satisfactory limit is hard to obtain. (2) A large condition number indicates a lack of knowledge 
in the inverse source determination process itself, so the predicted results will be questionable. (3) Any 
data errors or extraneous high frequency content present will be greatly exaggerated in the source estimates. 
These issues are clearly related to the associated Point Source Support Problem (PSSP). Regularization can 
be used to largely solve these conditioning problems in a natural way. 

To motivate what follows, consider a preliminary argument indicating a connection of mass dislocation 
and energy. Toward that end, consider the rather specialized situation where a very detailed reference model 
exists that fully represents the part of the mass density that is locally homogenous so that all that is left 
to predict are local density irregularities. Suppose, further, that this reference model has been subtracted 
from W and that a small enough uniform grid spacing is used so that one can directly associate a given 
point mass value with local density irregularities, which can be reinterpreted as a small block of matter. 
Two different physical scales will be considered, where the finer one is associated with this uniform grid of 
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point mass locations. (The point here is that generally one should distinguish between the fixed framework 
of point masses that are used for estimation purposes, which are generally assumed to be at fixed locations 
and the mass distributions that they are assumed to model, which may well shift.) For location k, let ruk 
be the point mass value in question and let Mk be the total mass of this block associated with the given 
potential W . If there are no local stresses in the medium, then the larger physical block of matter that 
Mk is part of is in its normal configuration and the subtracted reference field accounts for all of the local 
density pk so that mk = 0. Alternatively, consider what happens when the larger block of matter that nik 
is part of is subjected to compression along one direction, say the vertical direction. Let a capitol letter K 
be associated with this larger physical block of matter that Mk is part of so that its mass is represented by 
Mk- Then, let Lk, Wk and Dk denote the length, width and depth of this larger block. Further assume 
that the distortion of this larger block {SDk) is small and that the block responds in an elastic (i.e., linear) 
manner to this force by a change in Dk only. Since the block mass is conserved, Mk is constant and thus 
SpK = Mk/{Lk Wk [Dk — SDk]) — Mk/{Lk Wk Dk) ~ SDk Pk, so that m^ oc SpK oc SDk- Finally since 
the potential energy associated with elastic forces is proportional to {SDk)'^ it is clear that the energy of 
this internal dislocation is proportional to rn^ and thus the total energy for all the dislocations caused by all 
the various stresses can be written as X]fe=i '^k "m-k, where the afe's are constants of proportionality. If all of 
the blocks can be treated consistently, then this energy can be written as a X]fc=i "^fe- Adding this internal 
configuration energy to the external field energy form yields 

$= ||y-T4^||?, +an 

2 

as a more accurate replacement for $ = IIF — W^ll where 

f2 = ill := /^ iT^k ■ 

This is a standard quadratic regularization form that is invariably introduced solely on the grounds that it 
reduces the condition number. In particular, notice that using Cli effectively adds a diagonal term to Tk^k' 
and this clearly reduces the condition number of the linear equation set (which is especially obvious when 
normalized basis functions are used). 

_ 2 

The general philosophy underpinning the use of fi here is easily stated. When minimization of ||y — iy|| 

fails to specify a unique density estimate, the addition of fi will select those densities that have the lowest 
internal energy configuration, all other things being equal. (It is generally accepted that the action of a 
regularization form like fii produces a unique fit.) Finally before heading on it is worth discussing the 
implementation of normalized basis functions here. The action of a regularizing function is generally ignored 
when the normalization conditions are implemented so that it is fixed regardless of regularization. The 
regularization functions involving mass (like Oi), however, should be defined in terms of rfik rather than 
TOfc, but this complication is not considered in the text. (Further it worth noting that it will be assumed 
in the sequel that normalized basis functions are used; however, for convenience this normalization process 
is generally carried out without considering fi, and then the effects of this term arc added in just prior to 
computation of the linear equation set solution.) 

While the reduction of condition number associated with the use of f2i is desirable, there are three 
additional complications to consider here. First, in the above an elastic material medium that was surrounded 
by a like medium on all sides was considered — does a similar argument hold when these conditions do not 
hold? As an alternative example consider surface volume elements made up of a noncompressible fluid. As 
a realistic concrete example consider the ocean surface where a standard normal reference ellipsoid model 
has been subtracted off. Consider the following three well known facts [5]: (1) The sea surface is a surface of 
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constant potential. (2) The difference in altitude between this surface displacement and the normal ellipsoid 
is called geoid height (usually denoted H) . (3) The geoid height is proportional to the difference between the 
potential at the point in question on the surface and on the normal ellipsoid itself. Notice that while the 
density and sides of a surface volume element are fixed, the height varies and this change in height (the geoid 
height) leads directly to a change in mass that is proportional to the change in geoid height. (Alternatively 
for the point mass fitting algorithm the point mass values tend to be proportional to the potential difference 
that is fit, which is proportional in turn to the geoid height [in performing a point mass surface fit here 
for the region Qq, one generally chooses an appropriate Rq < Re so that the point masses do not tail-off 
away from Xj^ either too fast or to slow].) Thus here the energy is proportional to \mk\ rather than ?Ti^, 
so perhaps a better form of regularization would be to take aX]fc=i '^fc + f^'l2k=i 1"^^! ^^ ^ regularization 
factor. Perhaps even a factor proportional to X]fc=i l"^fe|^) where 1 < /i < 2, should be considered. Here the 
inclusion of a regularization term with 1 ^ /x < 2 leads to a nonlinear equation set for the m^'s, which is 
clearly inconvenient. Furthermore, since minimization of the form J7i tends to minimize X]fe=i I'^fcl'^ ^^ well, 
only the form Jli will be considered in the sequel (but this is clearly one of many numerous open issues). To 
summarize, when a linear restorative force is present (as one might expect for material stresses and strains) 
there is a well known quadratic dependence of (potential) energy, but when the displacement mechanism is 
directly related to the action of gravity on fluid surfaces, the potential energy tends to go like the well known 
mg H factor encountered in elementary physics books. 

The final two of the three objections to fJi can be stated briefly. The second one is that Cli tends to 
minimize the overall point mass values and if the subtracted reference model is not as detailed as required this 
will lead directly to systematic errors in the estimates. Thirdly and perhaps more importantly: Does even 
a first order approximation to the required reference model exist for cases of interest? These last two issues 
lead to a consideration of other regularization forms. There is a clear hint in the analysis performed above 
that leads to an improvement. In particular, since the stress arising on block k probably originated from a 
neighboring block in contact with it, a better model to consider is perhaps a regular grid of coupled blocks 
that can be viewed as a three dimensional assembly of masses that are coupled by springs in the vertical and 
horizonal directions. When such a system is in homogenous static equilibrium, the distortional energy is zero, 
but when each block is either compressed or stretched the resulting total energy will increase. Thus when one 
mass is displaced due to pressure from an adjacent mass, not only will the density of that particular block 
be increased, but also the energy of the block that is directly coupled to it. The configuration energy of such 
a coupled pair (just as for a coupled string configuration) can thus be represented by bJk^k'i'mk — rnk'Y ■ By 
minimizing this coupling energy, a smooth density profile results and the effect of the mass reduction effect 
is not as pronounced as it is with the straightforward regularization term fii. For ease of implementation 
here a quadratic form is desirable and it is also necessary to try to enforce (or at least strongly encourage) 
TOfc « TOfc/ for \X'j^ — X'j^i\ « 0. A regularization form that fulfills these requirements is 



Nk Nk 



Vl = Q.2:=2_^2^uJk,k'{'rnk-mk') (20) 

fc=i fe'=i 

where LOk,k' produces mass correlation effects. Thus, in general ojfc.fe' — ^k,k' — ^{dk,k') with dk^k' '■= 
\X'f. — X\.,\. Here, in particular, one suitable choice might be uj{d) — 1 if d is less than -s/S times the 
(three-dimensional) grid spacing and uj{d) = otherwise so that only the closest neighbors are correlated. 
This choice of regularization function reduces the condition number and introduces a uniformity into the vak 
values without reducing the overall mass values. This general form of regularization is desirable for many 
unrelated applications as well. With the right choice of ld and a, one can clearly negotiate very small grid 
spacing. 

Before considering several generally desirable refinements to this regularization process, notice that if 
one assumes that the variations in point mass values is proportional to the underlying local mass deficits 
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or excesses, and that these excesses and deficits are, in turn, proportional to an energy shift due to a local 
dislocation in the underlying material medium (which, by a standard Taylor's series argument, represents 
the displacement from what would otherwise have been a local energy minimum), then £72 is proportional 
to energy dislocation of the underlying material medium. This point will be elaborated on below. Mathe- 
matically (mk — ruk')^ = rrvl + ?n^, — 2mk'mk', so CI2 contains quadratic terms that effectively add to the 
diagonal of T^. k' as well a bilinear terms, which effectively subtract from the larger off-diagonal elements 
of Tk.k' (when |X^ — X^, | w Tk.k' ~ !)• Thus Ct2 actually has a stronger regularizing effect than Hi (all 
other things being equal). 

The foregoing regularization description might best be alternatively encapsulated in terms of a standard 
Taylor's series argument. Since $ = <i>(m) (where m is the vector of mass values), and in particular 
D, = n(m), it is natural to consider the energy variation in terms of m: 

g2 ^ 



dmkdmk' 



where Vm has components given by ■§■ — . In the simplest standard context that this Taylor series argument 
is used by physicists, an energy minimum, £, is sought and the displacement, x, is the variable of interest. 
Since it is argued that a physical minimum is present, linear terms cannot be present so the form 

£.{x) = £o+x^£" + ■■■ 

results (which ignores the possibility of physical terms of the form \x\). Since a minimum of <& is sought, the 
constant term for ^(m) can be ignored, and if it is assumed that a suitable reference has been subtracted 
off, m can be identified with 5xn; nevertheless, the linear terms obviously can not be ignored here. There 
are three reasons for this. First, fi is part of $ and linear terms definitely cannot be ignored in the rest of 
<& [see ([6]), for example] since the linear terms might cancel out in some fashion. Second, the regularizing 
functions discussed below (i.e., X][-^o — "^fc]^) have linear factors as well as bi-linear terms (X]["^fc — m-fe']^); 
moreover, these regularization forms have been shown to posses physical relevance. Third, as noted above, 
regularization factors proportional to ^ \mk\ have a reasonable physical basis. At any rate, the physical 
significance of the first few Taylor series terms should be apparent. Finally note that, with respect to 
incompressible fluids and/or stratified media, the above regularization analysis is incomplete at best. Here 
the most relevant factor is the shape of the media boundary surface separating one density layer type form 
another (consider, for example, the ocean floor). (Hopefully some sort of future analysis undertaken by others 
will demonstrate a more reflned understanding of the higher order aspects and of appropriate regularization 
functions, in general.) A means of meshing regularization and surface boundary information will thus be 
considered next. 

Suppose there are various regions, or density layers, which are distinct, but that each such region, or 
layer, tends to be homogenous. This situation can easily be modeled by using a proper choice of uJk^k' in 
0,2- Thus let {'R-.j}j;ii for J = 1, 2, 3, . . . , Nj be a suitable partition of fl'^ into subregions: 

Nj 

\Jnj = n', and n.j n nj> =0 for j 7^ j' 

where J' = 1, 2, 3, ... , Nj. Then let 

i^k.k' — ^{dk,k') if X'l^ and X^., e TZj, and otherwise let ujk,k' = 0. (21) 

The resulting regularization approach characterized by (|21|) tends to produce independent homogenous den- 
sities for each of the separate regions. 
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Alternatively, suppose that rough density profile information is available from seismology (or some other 
means) and this profile is specified by p{X'), then this profile can be discretized: p(X') =^ rhk- Then a 
regularization term of the form 

rt ^0,3 := ^^{mk - rhkf 

k=l 

may be appropriate. One special case is rhk = Mj for X' e TZj, where Mj is a constant. This general 
type of regularization is clearly appropriate when density variations are highly depth dependent and overall 
density averages are known for certain depths, such as the core or deep mantle. 

Structured Point Source Technique (SPST) Basis Functions 

This subsection first briefly recaps some of the points made in Section (ii), where the general ideal of a 
SPST basis function was introduced, and then gives a mathematical description. 

Thus, as previously noted, a common gravitational source density prediction strategy is to divide up 
the source region into a collection of regular homogenous bodies — with the most common example being to 
divide up fl'g into arrangement of (nonintersecting) regular parallelepipeds that come close to covering the 
entire region of interest. This has an intrinsic regularizing effect. A very flexible way to implement such a 
scheme in the present formalism involves approximating the field of each such subregion by a regular grid 
of closely spaced point masses. When a collection of point masses share a single common constant mass 
value (or have fixed relative mass values) and thus have only one undetermined source term, the resulting 
structure will be labeled a structured point source. In order to better understand this from a regularization 
perspective, it is useful to compare this to some of the regularization schemes just considered. Clearly Q2 
generally has the effect of forcing ruk ~ m-fe' for nearby masses ruk and mk' ; moreover, when it is implemented 
according to the partitioned region regularization scheme, as characterized by using (j2ip . similar end results 
to a SPST basis function fit might well be expected (although one might reasonably argue that SPST basis 
functions have a stronger regularization effect). Next consider the regularization effects of ^3 versus SPST 
basis functions. Again one might expect similar end results for most implementations, but clearly, when all 
else is equal, the automatic constraint implicit in the structured point source technique will have a stronger 
regularizing effect. 

As above, let {'R'j}j;Li be a partition of il(. into subregions: 

Nj 



y 7^J = n'^ and Uj n 7^J' = for J :^ J' 



J=l 

where J = 1, 2, 3, ... , Nj and J' = 1, 2, 3, ... , Nj. Then the basic SPST idea is to hold the mass fixed 
over each subregion: m^, := Mj for all k such that X'^, g TZj. One advantage of this approach over the 
regularization approach described by (|20p and (j2ip is that the resulting linear equation set has dimension 
Nj X Nj rather than Nk x Nk- In detail, let j be a local index for each of the TZj and let n{J) be the number 
of (uniform) point sources in TZj, then X)j=i "-('^) — -^fc ^^^ 1 = J = n{J). Further let K{J, j) denote a 
reordering of the index k such that for all 1 ^ j ^ "i^iJ) and 1 ^ J ^ Nj 

X'k{J. ]) e -Rj . 
Then the resulting point mass potential field can be written 

Nj n{J) 

Vn, = y^ Mj\^ — i^ 

J=l , = 1 ' ^(■^■^■)' 
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and the resulting DIDACKS SPST linear equation set for the Mj's is 

Nj 



J2 ^j. J' ^''' = ^j 



j'=i 

The required expressions for the half-space energy norm are 

^(J) 



^^ = iE^(^^w.)) 



and 
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U,J' 



n{J) n(J') 

7 / . / . rp —^, 1 • Here Pk{j,j) = {^'k{j,,y v'k(j,j). ~^'k(j,j)) 



T 



— ; , , ^K(.j,3)-^K{j'.r) 

J = 1 J = 1 y J I 



Observe that the above description of SPST basis functions can be easily generalized along the lines 
discussed in Section (ii) by introducing fixed relative scaling factors. The idea, of course, is to fix the relative 
ratios of the various component point masses and then determine the overall scale of the configuration by 
the fitting process. Hence, for a set of point masses nik, using the given rations C/t, set nik = CkMj, for 
Xj. € TZj, where only the Mj need to be determined and these values can be determined much the same as 
before. This idea has a number of potentially useful applications along the lines of the examples mentioned at 
the end of Section (ii). Some implementation points, however, may not be completely transparent. Suppose, 
for example, that the boundary of the regions TZys are not known and that one wants to estimate these 
boundaries by a means of DIDACKS NLLSQ scheme. It is easy enough to set up the problem and get an 
appropriate cost function to use, but there is an underlying fixed grid for the point sources locations, so that 
when the surface moves a source point will abruptly switch from one region to another. A NLLSQ algorithm 
generally requires good partial derivative information to work well (although there are specialized discrete 
optimization approaches) so there is a problem. Thus the problem is that the point mass grid spacing will 
generally be such that a perturbation of the boundary may not change the cost function, or it may change 
abruptly all at once, so that the resulting partial derivative information will not be acceptable for NLLSQ 
purposes. The easiest solution to this problem is to use the above idea involving point mass ratios in order 
to implement "soft boundaries." This is easily accomplished as follows. Although the details may not be 
of interest to many readers, they are included here since they show some of the power and flexibility of the 
SPST basis function concept. 

For concreteness suppose that the half-plane region Oi is of interest and that the boundaries of TZj are 
generally parallel to the x', y' and z' axes. In this case it is especially easy to parameterize the boundary 
surfaces of TZj. Observe that, because of shared boundaries, when TZj is surrounded on all sides one can, for 
accounting purposes, assume that three of the six sides are associated with TZj, so that there are nominally 
six values of {rji} associated with each TZj. 

The overall point is to keep the idea of a well delineated boundary and to keep a fixed SPST basis 
function for each TZj, but to define the basis functions in such a way as to make a smooth transition in 
mass profile across the boundary, without leaving gaps — that is, if the given potential field is constant and 
uniform over a prediction region, then the resulting mass estimates should be too for this case. For purposes 
of this definition, let each region TZj and its adjacent neighbor TZji overlap by a certain distance D. The 
idea is that over this distance a gradual transition is made from the uniform mass of TZj to the uniform mass 
of TZji. For example, if TZj and its adjacent neighbor TZji share a common face along the x— axis then the 
X coordinates for the transition might be labeled X!p ^ x' ^ X!j, + D where X!p = constant -I- rji for some 
i (and the constant here is taken so that the midpoint of this interval corresponds to the boundary of TZj 
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and its adjacent neighbor TZj'). Then if Xk is a point in this transition zone, when X'^, is treated as a point 
in TZj one should set the corresponding value of the mass for this region to /i(xjj — X'rp)Mj where h{x) is 
a fifty percent cosine taper or h{x) — x/ D. [A fifty percent cosine taper for the interval [o, 27r] is simply 
(1 + cosa;)/2.] When the same point is considered as a point in TZji then the corresponding mass value is 
set to [1 — h{x'j^ ~ X'rp)]M'j so that there is a gradual transition from TZj to TZ'j. When the DIDACKS cost 
function with these SPST basis functions is set up, it will have the form '^{M.j'q), where rj is the vector of 
the entire set of surface parameters. Minimizing this cost function for both Ad,j and the components of 77 
requires a NLLSQ algorithm for its solution. 

Next consider the second application of the above mass ratio concept. Although the following idea is 
general, for the sake of simplicity consider the case where there is only one region of interest so that TZi := fJ^ 
and 



VnAX) =^ 



TOfc 



\x-x',\ 



The central idea here is that instead of using mk directly as a fitting parameter, new set of mass fitting 
parameters is introduced: 

Nc 
ruk = nikiCi) = ^Ci*,(Xfe) 

i=l 

where the 4'i's are a set of Nc suitable basis functions. The minimization of the resulting cost function 
$ yields a linear set of Nc equations for the C^'s. Generally here one might place the X'j, on a (tight) 
three dimensional uniform grid; however, it is also possible to arrange the X^'s on a surface grid or along 
a line array. When the X^'s are arranged on a line a good choice for the ^i's might be a set of Fourier 
series basis functions expressed as functions of path length along the line. In these sort of approaches 
regularization should still be applied as needed. Also observe that one can easily extend this idea to separate 
basis expansions over each of the subregions TZj, so no real loss of generality resulted from considering the 
special case TZi :— il'^. Finally, the actual form of the linear equation sets that result for these discrete 
parameterized fits are easily written down and the implementation details are straightforward. Moreover, 
the resulting equations are quite similar to those that result from parameterized continuous distributions, 
which we now turn to (where integrals just basically replace sums). 

Continuous Source Estimation 

It is a small step from the discrete parameterized fits just considered to the consideration of full continuous 
distributions. Continuous distributions require numerical integration, so they are more difficult to implement. 
To streamline the presentation only the ili case will be considered here (the J7o case follows in a like fashion). 
Here 

ViX) . 




where pv is a parameterized continuous distribution, which can be taken to have the following form 



Nc 



pM')^TCnMX') 
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Here again, the ^„'s are a set of suitable basis functions. The resuhing hnear equation set for the C„'s is 

Nc 



where 



"'" ^ J J J J J J Vi^' - ^")' + iy' - v"? + (2' + z"? 



with X' and X" G fl' Likewise 



^" = \jjl W{x', y', -z')^{x', y', ~z') d' X' 
n 
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